; docformat = 'rst'
;
;+
;
; :Purpose:
;   Process MOZART output and calculate EUM sensitivites for period before start year
;   
; :Inputs:
;   sim_name: (string) simulation name
;   
; :Requires:
;   Directory of MOZART outputs to be given in {sim_name}.params file. 
;   Output files must contain MOZART variables corresponding to the IC_NAME parameter
;
; :Outputs:
;   Generates {input_directory}/{sim_name}/EUM/sensitivity_spinup.inp.sav file
;   
; :History:
;   Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
pro cels_mozart_sensitivity_spinup, sim_name

  compile_opt idl2
  on_error, 2
  
  print, "CELS_MOZART_SENSITIVITY_SPINUP"

  Lev_Index=-3  ;Level retrieve from MOZART
  StartY=cels_get_parameter(sim_name, 'STARTY')
  StartY_IC=cels_get_parameter(sim_name, 'STARTY_IC')
  Stations=cels_get_parameter(sim_name, 'STATIONS')
  Short_name=cels_get_parameter(sim_name, 'POLLUTANT_SHORT')
  Stations_all=[stations, cels_get_parameter(sim_name, 'BG_STATIONS')]

  IC_names=cels_get_parameter(sim_name, 'IC_NAME')
  EUM_OUTPUT_DIRECTORY=cels_get_parameter(sim_name, 'EUM_SPINUP_DIRECTORY')

  nStations=n_elements(stations)
  nStations_all=n_elements(stations_all)

  ;Get simulation info
  restore, cels_filestr(/Input, sim_name + '/measurements/measurements_spinup.sav')
  restore, cels_filestr(/Input, sim_name + '/state.sav')
  restore, cels_filestr(/Input, sim_name + '/regions.sav')
  restore, cels_filestr(/Input, sim_name + '/emissions/emissions_totals.inp.sav')

  ;Variable names to extract from MOZART
  EUM_ID_spinup=IC_NAMES

  ;Check station ordering
  undefine, mf
  if total(mf_station_ID ne stations_all) then begin
    message, "Measurement codes don't match"
  endif
  
  EUM_SPINUP_STATIONS=STATIONS_ALL

  files=file_search(EUM_OUTPUT_DIRECTORY, '*h0*')
  nFiles=n_elements(files)

  ;Get output file grid
  lon=ncget('lon', filename=files[0])
  lat=ncget('lat', filename=files[0])
  if total(lon ne lon_eum) gt 0. then begin
    Print, "CELS_MOZART_SENSITIVITY_SPINUP: WARNING, MOZART spinup grid doesn't match grid.inp.sav"
  endif
  lonsize=n_elements(lon)
  latsize=n_elements(lat)

  ;Measurement locations
  slonI=intarr(nstations_all)
  slatI=intarr(nstations_all)
  for si=0, nStations_all-1 do begin
    dummy=min(lon - mf_lon[si], LonI, /abs)
    dummy=min(lat - mf_lat[si], LatI, /abs)
    slonI[si]=lonI
    slatI[si]=latI
  endfor


  ;Declare output arrays
  timesteps_in_eum=12L*Long(StartY - StartY_IC) ;Assumes MONTHLY model output
  eum_spinup_time=dblarr(timesteps_in_eum)
  H_Spinup=fltarr(N_IC, nStations_all, timesteps_in_eum)

  timesteps_cumulative=0L

  ;Get model output
  for fi=0, nfiles-1 do begin
    time_file=mozdate(ncget('date', filename=files[fi]))
    wh_time=where(time_file ge julday(1, 1, StartY_IC, 0) and time_file le julday(1, 1, StartY, 0), timesteps_in_file)

    if timesteps_in_file ne 0 then begin
      eum_spinup_time[timesteps_cumulative:timesteps_cumulative+timesteps_in_file-1]=time_file[wh_time]

      For xi=0, N_IC-1 do begin
        var_Pert=ncget(short_name + '_' + EUM_ID_spinup[xi] + '_VMR_avrg', filename=files[fi])   ;MLR REMOVE UNDERSCORE!
        var_Pert=reform(var_Pert[*, *, lev_index, wh_time])
  
        ;Assumes no reference run (i.e. perturbation relative to zero)
        H_global=var_Pert/Pulse_magnitude_IC[xi]
        for si=0, nStations_all-1 do begin
          H_spinup[xi, si, timesteps_cumulative:timesteps_cumulative+timesteps_in_file-1]=$
            H_global[slonI[si], slatI[si], *]
        endfor
      endfor

      timesteps_cumulative+=timesteps_in_file
    endif else begin
      if fi eq 0 then begin
        print, "CELS_MOZART_SENSITIVITY_SPINUP: WARNING, first file outside of time range"
      endif
    endelse;time in file

  endfor
  
  H_spinup=h_spinup*1.e12
  
  caldat, eum_spinup_time-15.d, month, dummy, year
  eum_spinup_time=julday(month, 15, year, 0)  ;Assymes monthly

  save, filename=cels_filestr(/Input, sim_name + '/EUM/sensitivity_spinup.inp.sav'), $
    eum_spinup_time, H_spinup, EUM_ID_spinup, EUM_SPINUP_STATIONS
 
End